setwd("/mnt/md1200/6/yjp/5hmc_analysis_hg19_new/20210102")

filea=read.csv("../20201207/all.FDR.sig.at.least.one.add.direction.same.diff.csv",header=T)
filea$id=paste(filea$Chr,filea$Start,sep = ":")
filea1=filea[filea$FDR.sig>1,]

file=read.csv("../20201207/at.least.one.AShM.in.DC.add.BF.beta0.add.CCHC.csv",header=T)
file$id=paste(file$Chr,file$Start,sep=":")
file1=file[file$pattern.not.rm.dupl.num.DC>1,]
file2=file1[file1$BF_in_DC>1,]
file3=file1[file1$BF_in_DC>10,]

#####motifbreakR
library(motifbreakR)
library(BSgenome.Hsapiens.UCSC.hg19)
library(SNPlocs.Hsapiens.dbSNP142.GRCh37)
library(SNPlocs.Hsapiens.dbSNP144.GRCh37)
library(BiocParallel)
												###pctRef和pctAlt显示序列中基序的分数，作为基序在理想序列上可以达到的最佳分数的百分比。
												#参考：http://bioconductor.org/packages/devel/bioc/vignettes/motifbreakR/inst/doc/motifbreakR-vignette.html
snps=as.character(filea$avsnp150)
snps=snps[!snps=="."]

search.genome=BSgenome.Hsapiens.UCSC.hg19
dbSNP=SNPlocs.Hsapiens.dbSNP142.GRCh37
rsid=snps

change.to.search.genome <- function(granges.object, search.genome) {
  if (Reduce("&", !is.na(genome(granges.object)))) {
    if (identical(genome(granges.object), genome(search.genome))) {
      return(granges.object)
    }
  }
  if(isTRUE(all.equal(seqlevels(granges.object), seqlevels(search.genome)))) {
    seqinfo(granges.object) <- seqinfo(search.genome)
  } else {
    if(seqlevelsStyle(granges.object) != seqlevelsStyle(search.genome)) {
      seqlevelsStyle(granges.object) <- seqlevelsStyle(search.genome)
    }
    normal.xome <- seqlevels(granges.object)[(regexpr("_", seqlevels(granges.object)) < 0)]
    #xome.value <- str_extract(normal.xome, "[0-9]|1[0-9]|2[0-9]|3[0-9]|4[0-9]|5[0-9]|6[0-9]|7[0-9]|8[0-9]|9[0-9]|X|Y|M")
    positions <- unlist(sapply(paste0(normal.xome, "$"), grep, seqnames(seqinfo(search.genome))))
    new2oldmap <- rep(NA, length(seqinfo(search.genome)))
    new2oldmap[positions] <- 1:length(positions)
    seqinfo(granges.object, new2old = new2oldmap) <- seqinfo(search.genome)
  }
  return(granges.object)
}

rsid.grange <- as(snpsById(dbSNP, rsid, ifnotfound = "warning"),"GRanges")
rsid.grange <- change.to.search.genome(rsid.grange, search.genome)
rsid.grange <- GRanges(rsid.grange)
rsid.refseq <- getSeq(search.genome, rsid.grange)
rsid.grange$UCSC.reference <- as.character(rsid.refseq)
rsid.grange <- split(rsid.grange, rsid.grange$RefSNP_id)

determine.allele.from.ambiguous <- function(ambiguous.allele, known.allele) {
  neucleotide.ambiguity.code <- list(Y = c("C", "T"), R = c("A", "G"), W = c("A", "T"),
                                     S = c("G", "C"), K = c("T", "G"), M = c("C", "A"),
                                     D = c("A", "G", "T"), V = c("A", "C", "G"),
                                     H = c("A", "C", "T"), B = c("C", "G", "T"),
                                     N = c("A", "C", "G", "T"))
  specnac <- neucleotide.ambiguity.code[[ambiguous.allele]]
  unknown.allele <- specnac[-grep(known.allele, specnac)]
  return(unknown.allele)
}
snps=rsid.grange
rt.alt.allele=c()
rsid2=c()
for (i in rsid){
snp=snps[[i]]
if (!is.null(snp)){
alt.allele <- determine.allele.from.ambiguous(snp$alleles_as_ambig, snp$UCSC.reference)
if (length(alt.allele)>=1){
rt.alt.allele=c(rt.alt.allele,alt.allele)
rsid2=c(rsid2,i)
}
}
}
length(rsid2)

snps=rsid2

snps.mb <- snps.from.rsid(snps,dbSNP = SNPlocs.Hsapiens.dbSNP142.GRCh37,search.genome = BSgenome.Hsapiens.UCSC.hg19)
pwmlist=query (query (MotifDb, 'hsapiens'), 'SELEX')
rt=motifbreakR(snps.mb, pwmList=pwmlist, threshold = 1e-4, filterp = TRUE,method = "log", show.neutral = FALSE, verbose = FALSE,bkg = c(A = 0.25, C = 0.25, G = 0.25, T = 0.25), BPPARAM = bpparam())
a=data.frame(rt)
a$rsid=names(rt)
write.table(a,"117K.ASH.motif.predict.txt",quote=F,row.names=F,sep="\t")

